LSTM AutoEncoder¶
気象情報をセンサデータと見立てて異常検知してみる。入力は気象庁の過去の天気情報で入手。1980年から2020年(06月まで)の毎日の天気情報。
In [1]:
import pandas as pd
import numpy as np
import os
from tqdm import tqdm
import tensorflow as tf
import pandas as pd
pd.options.mode.chained_assignment = None
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
import plotly.express as px
import plotly.graph_objects as go
%matplotlib inline
sns.set(style='whitegrid',palette ='muted')
rcParams['figure.figsize']=14, 8
np.random.seed(1)
tf.random.set_seed(1)
print("TensorFlow Version: ",tf.__version__)
データの準備¶
予め超かんたんなクレンジングを実行済み。input配下にCSVが入れてある状況。
サンプルとして一番初めの月だけ見てみる。その後全期間のデータをマージ、風向きデータのOneHotが面倒なので可視化後に除外。数値データだけ残して正規化する。
※最高気温、最低気温、平均気温など相関の強いカラムが多いので本来はカラムを絞ったり次元削減するべきだが、とりあえずのモデルなのでそこまではこだわらない。
In [2]:
INPUT_DIR = 'input/'
In [3]:
csv_files = os.listdir(INPUT_DIR)
csv_files.sort()
In [4]:
# とりあえず初めの月だけ表示してみる
input_path = os.path.join(INPUT_DIR, csv_files[0])
df = pd.read_csv(input_path)
df
Out[4]:
In [5]:
# 全データを結合
df = pd.DataFrame(columns=df.columns)
for csv_file in tqdm(csv_files):
input_path = os.path.join(INPUT_DIR, csv_file)
tmp_df = pd.read_csv(input_path, parse_dates=['date'])
df = pd.concat([df, tmp_df])
df = df.fillna('').reset_index(drop=True)
In [6]:
# 風向きはいらない
direction_cols = ['wind_velocity_max_direction', 'maximum_instantaneous_wind_speed_direction', 'most_common_wind_direction']
train_cols = [col for col in df.columns.tolist() if not (col in direction_cols)]
df = df[train_cols]
In [7]:
train_size = int(len(df)*0.85)
test_size = len(df)-train_size
train,test = df.iloc[0:train_size],df.iloc[train_size:len(df)]
display(train.tail(5))
display(test.head(5))
print(train.shape,test.shape)
feature_cols = train.drop('date', axis=1).columns.tolist()
可視化¶
超大雑把に可視化しておく。残念ながら積雪や降雪の記録は取れていないらしい。三浦地方は大雪は少ないから異常がわかりやすいと思ったのに。。。
pandas-profilingがものすごくリッチに成長していてビックリ。先述の通り相関が高いものを警告したりしてくれてる。
In [8]:
for col in feature_cols:
fig = go.Figure()
fig.add_trace(go.Scatter(x=df['date'],y=df[col],mode='lines',name=col))
fig.update_layout(showlegend=True)
fig.show()
In [9]:
import pandas_profiling as pdp
pdp.ProfileReport(df)
Out[9]:
正規化¶
標準化、正規化は悩みどころだけど、あとで各カラムの予測値と実測値の差分を足し合わせる関係で正規化にしておく
In [10]:
from sklearn.preprocessing import StandardScaler, MinMaxScaler
X_train = train[feature_cols]
X_test = test[feature_cols]
scaler = MinMaxScaler()
scaler = scaler.fit(X_train)
X_train = scaler.fit_transform(X_train)
X_test = scaler.fit_transform(X_test)
X_train = X_train.reshape(X_train.shape[0], 1, X_train.shape[1])
X_test = X_test.reshape(X_test.shape[0], 1, X_test.shape[1])
モデル生成¶
あまり深くする必要もないのだが、せっかくなので多少は積んでおく。推論が遅くて使えなかったら考える。
In [11]:
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import Input,LSTM,Dense,Dropout,RepeatVector,TimeDistributed,BatchNormalization
from tensorflow.keras import regularizers
def build_lstm_ae(X_train):
inputs = Input(shape=(X_train.shape[1], X_train.shape[2]))
x = LSTM(256, activation='relu', return_sequences=True,
kernel_regularizer=regularizers.l2(0.00))(inputs)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(128, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(64, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(32, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(16, activation='relu', return_sequences=False)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = RepeatVector(X_train.shape[1])(x)
x = LSTM(16, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(32, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(64, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(128, activation='relu', return_sequences=True)(x)
# x = Dropout(0.1)(x)
x = BatchNormalization()(x)
x = LSTM(256, activation='relu', return_sequences=True)(x)
output = TimeDistributed(Dense(X_train.shape[2]))(x)
model = Model(inputs=inputs, outputs=output)
return model
In [12]:
model = build_lstm_ae(X_train)
model.compile(loss='mae',optimizer='adam',metrics=['accuracy'])
model.summary()
In [13]:
# AuroEncoderなのでラベルは不要。特徴量とラベルがおんなじって感じ。仮にノイズ除去モデルを作りたいなら第一引数をノイズを付与したX_trainにする
es = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5, mode='min')
# sb = tf.keras.callbacks.ModelCheckpoint("tf_model.hdf5", monitor="val_loss", verbose=1, save_best_only=True, save_weights_only=False)
history = model.fit(
X_train, X_train, epochs=100,
batch_size=256, validation_split=0.1, callbacks=[es], shuffle=False
)
In [14]:
model.save('tf_model.hdf5')
In [15]:
history.history.keys()
# fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False)
# axes[0].plot(history.history['loss'],label='Training Loss')
# axes[0].plot(history.history['val_loss'],label='Validation Loss')
# axes[1].plot(history.history['accuracy'],label='Training ACC')
# axes[1].plot(history.history['val_accuracy'],label='Validation ACC')
# plt.legend()
Out[15]:
In [16]:
plt.plot(history.history['loss'],label='Training Loss')
plt.plot(history.history['val_loss'],label='Validation Loss')
plt.legend()
Out[16]:
In [17]:
plt.plot(history.history['accuracy'],label='Training ACC')
plt.plot(history.history['val_accuracy'],label='Validation ACC')
plt.legend()
Out[17]:
In [18]:
model.evaluate(X_test,X_test)
Out[18]:
異常検知¶
AutoEncoderなので予測値は実際の気象情報を変換して復元した値になる。復元した値と実際の値の差分を取って、差分が大きいところを異常とする。
気象情報なので多少は異常があると仮定しているが、仮に本物のセンサデータで異常発生時のデータが無い場合は差分の平均との差分が+3σより大きかったら異常とか、観測データの差分の最大値+αを超えたら異常って扱いにするのかな。
差分はMAEとMSEで見ているが、正規化しているから差分は1未満になるのでMAEよりMSEのほうが顕著に出ているように見える。
BatchNorm入れたらかなり再現力が高くなって分布が見づらくなった。。。分布図をLogスケールで見る方法って無いのかな。。。
In [19]:
# plot the loss distribution of the training set
X_pred = model.predict(X_train)
X_pred = X_pred.reshape(X_pred.shape[0], X_pred.shape[2])
X_pred = pd.DataFrame(X_pred, columns=feature_cols)
X_pred.index = train.index
scored = pd.DataFrame(index=train.index)
Xtrain = X_train.reshape(X_train.shape[0], X_train.shape[2])
scored['Loss_mae'] = np.mean(np.abs(X_pred-Xtrain), axis = 1)
scored['Loss_mse'] = np.mean(np.power(X_pred-Xtrain, 2), axis=1)
plt.figure(figsize=(16,9), dpi=80)
plt.title('Loss Distribution', fontsize=16)
sns.distplot(scored['Loss_mae'], bins=100, kde=True, color='red', label='mae')
sns.distplot(scored['Loss_mse'], bins=100, kde=True, color='blue', label='mse')
plt.xlim([0.0,max(scored['Loss_mae'].max(), scored['Loss_mse'].max())])
Out[19]:
In [20]:
scored
Out[20]:
In [21]:
X_pred = model.predict(X_test)
X_pred = X_pred.reshape(X_pred.shape[0], X_pred.shape[2])
X_pred = pd.DataFrame(X_pred, columns=feature_cols)
X_pred.index = test.index
scored = pd.DataFrame(index=test.index)
Xtest = X_test.reshape(X_test.shape[0], X_test.shape[2])
scored['mae'] = np.mean(np.abs(X_pred-Xtest), axis = 1)
scored['mse'] = np.mean(np.power(X_pred-Xtest, 2), axis = 1)
plt.figure(figsize=(16,9), dpi=80)
plt.title('Test Error Distribution', fontsize=16)
sns.distplot(scored['mae'], bins=100, kde=True, color='red', label='mae')
sns.distplot(scored['mse'], bins=100, kde=True, color='blue', label='mse')
plt.xlim([0.0,max(scored['mae'].max(), scored['mse'].max())])
Out[21]:
In [22]:
X_test_pred = model.predict(X_test)
X_test_pred = X_test_pred.reshape(X_test_pred.shape[0], X_test_pred.shape[2])
Xtest = X_test.reshape(X_test.shape[0], X_test.shape[2])
mae = np.mean(np.abs(X_test_pred-Xtest), axis = 1)
mse = np.mean(np.power(X_test_pred-Xtest, 2), axis = 1)
test_df = test.copy()
test_df['mse'] = mse
test_df['mae'] = mae
In [23]:
# 全体のthreshold*100 % は正常値として扱い、それ以外を異常とする。そのための閾値を計算する。
threshold = 0.998
mse_threshold = np.quantile(mse, threshold)
mae_threshold = np.quantile(mae, threshold)
print('MSE {0} threshold:{1}'.format(threshold, mse_threshold))
print('MAE {0} threshold:{1}'.format(threshold, mae_threshold))
In [24]:
display(test_df[test_df['mse'] > mse_threshold])
display(test_df[test_df['mae'] > mae_threshold])
In [25]:
test_df['mse_threshold'] = mse_threshold
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse'],mode='lines',name='mse'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse_threshold'],mode='lines',name='mse_threshold'))
fig.update_layout(showlegend=True)
fig.show()
test_df['mae_threshold'] = mae_threshold
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae'],mode='lines',name='mae'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae_threshold'],mode='lines',name='mae_threshold'))
fig.update_layout(showlegend=True)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse'],mode='lines',name='mse'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse_threshold'],mode='lines',name='mse_threshold'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae'],mode='lines',name='mae'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae_threshold'],mode='lines',name='mae_threshold'))
fig.update_layout(showlegend=True)
fig.show()
実際にはリアルタイムに流れてくるデータを定期的に刈り取って予測すると思われるので、テストデータを定期的に区切って予測してみる。特に意味は無いけど、30件までしか一度に予測しない前提で行う。
あれ、頂点の値がさっきと一緒？LSTMってそういうものだっけ。時系列を見るから古いデータを捨てて予測したら別の値が出るものだと思ってたけど。。。
In [26]:
X_test.shape
Out[26]:
In [27]:
model = tf.keras.models.load_model('tf_model.hdf5')
In [28]:
offset = 0
for i in range(2200):
fetch_size = np.random.randint(4,10)
offset = offset + fetch_size
X_test_tmp = X_test[0:offset][-30:]
X_test_pred_tmp = model.predict(X_test_tmp)
X_test_pred_tmp = X_test_pred_tmp.reshape(X_test_pred_tmp.shape[0], X_test_pred_tmp.shape[2])
Xtest_tmp = X_test_tmp.reshape(X_test_tmp.shape[0], X_test_tmp.shape[2])
mae = np.mean(np.abs(X_test_pred_tmp-Xtest_tmp), axis = 1)
mse = np.mean(np.power(X_test_pred_tmp-Xtest_tmp, 2), axis = 1)
test_df = test[0:offset][-30:].copy()
test_df['mse'] = mse
test_df['mse_threshold'] = mse_threshold
if test_df['mse'].max() > mse_threshold:
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse'],mode='lines',name='mse'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse_threshold'],mode='lines',name='mse_threshold'))
fig.update_layout(showlegend=True)
fig.show()
In [ ]:
In [ ]: